Introduction

This document shows how to calculate candidate Pillars and their constituent Indicators at Census Block Group, Tract, and Service Area/ Municipality level. Here, we load the data we pre-processed in the previous exercise

b <- read_sf("../data/boundary.geojson")
bg.l <- read_sf("../data/blockgroups_long.geojson")
bg.w <- read_sf("../data/blockgroups_wide.geojson")
tr.l <- read_sf("../data/tracts_long.geojson")
tr.w <- read_sf("../data/tracts_wide.geojson")
pl.w <- read_sf("../data/place_wide.geojson")
pl.l <- read_sf("../data/place_long.geojson")
lead <- read_sf("../data/parcels_lead_service_lines.geojson")
lead.c <- sf::st_centroid(lead)

calculate_volume <- function(hh_size = 4, gpcd = 50){
  vol = hh_size * gpcd * 30
  return(vol)  
}

calculate_bill <- function(volume = 6000){
  fixed_charges = 7.63 + 9.85 + 1.8
  vol_rate = 5.76 + 2.71 
  bill = fixed_charges+ (vol_rate*volume)/748.052
  return(bill)
}



gradeAffordability <-function(HBI,PPI){
 grade <- NA
 grade <- case_when(
    HBI < 7 & PPI < 20 ~ "Low Burden",
    HBI < 7 & PPI>=20 & PPI <35 ~ "Moderate-Low Burden",
    HBI < 7 & PPI >=35 ~ "Moderate-High Burden",
    HBI >= 7 & HBI < 10 & PPI >= 35 ~ "High Burden",
    HBI >= 7 & HBI < 10 & PPI >= 20 & PPI < 35 ~ "Moderate-High Burden",
    HBI >= 7 & HBI < 10 & PPI < 20 ~ "Moderate-Low Burden",
    HBI >= 10 & PPI >= 35 ~ "Very High Burden",
    HBI >= 10 & PPI >= 20 & PPI < 35 ~ "High Burden",
    HBI >= 10 & PPI < 20 ~ "Moderate-High Burden"
  )
  return(grade)
  
}


save.image("../cache/components_prep.RData")
load("../cache/components_prep.RData")

Pillar 1 | Access: Affordability and Service Quality

1.1 Affordability

1.1.1 HBI/PPI Matrix

We adopt the approach in the figure below, as elaborated in this report.

HBI/PPI Matrix

Overall (Census Place - level data)
#HBI Calculation

# Place
pl.w <- pl.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
pl.w <- pl.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
pl.w <- pl.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
pl.w <- pl.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
pl.w <- pl.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
pl.w <- pl.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
pl.w <- pl.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
pl.w <- pl.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))


pl.w$HBI_size_avg <- 100*pl.w$bill_size_avg/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_1 <- 100*pl.w$bill_size_1/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_2 <- 100*pl.w$bill_size_2/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_3 <- 100*pl.w$bill_size_3/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_4 <- 100*pl.w$bill_size_4/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_5 <- 100*pl.w$bill_size_5/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_6 <- 100*pl.w$bill_size_6/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_7 <- 100*pl.w$bill_size_7/(pl.w$hh_income_upper_limit_Q1E/12)


pl.w$PPI <- 100*(pl.w$inc_poverty_countE-pl.w$inc_poverty_count_gr2E)/pl.w$inc_poverty_countE
pl.w$aff_grade_avg <- gradeAffordability(pl.w$HBI_size_avg,pl.w$PPI)
pl.w$aff_grade_1 <- gradeAffordability(pl.w$HBI_size_1,pl.w$PPI)
pl.w$aff_grade_2 <- gradeAffordability(pl.w$HBI_size_2,pl.w$PPI)
pl.w$aff_grade_3 <- gradeAffordability(pl.w$HBI_size_3,pl.w$PPI)
pl.w$aff_grade_4 <- gradeAffordability(pl.w$HBI_size_4,pl.w$PPI)
pl.w$aff_grade_5 <- gradeAffordability(pl.w$HBI_size_5,pl.w$PPI)
pl.w$aff_grade_6 <- gradeAffordability(pl.w$HBI_size_6,pl.w$PPI)
pl.w$aff_grade_7 <- gradeAffordability(pl.w$HBI_size_7,pl.w$PPI)


pl.bill<- pl.w %>% st_drop_geometry() %>% select(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7) %>%
     pivot_longer(cols = c(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7),names_to="hh_size",names_prefix="bill_size_",values_to="bill")

pl.HBI<- pl.w %>% st_drop_geometry() %>% select(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7) %>%
     pivot_longer(cols = c(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7),names_to="hh_size",names_prefix="HBI_size_",values_to="HBI")

pl.aff<- pl.w %>% st_drop_geometry() %>% select(PPI,hh_income_upper_limit_Q1E,aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7) %>%
     pivot_longer(cols = c(aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7),names_to="hh_size",names_prefix="aff_grade_",values_to="aff_grade")

pl.aff <- pl.aff %>% left_join(pl.HBI,by="hh_size") %>% left_join(pl.bill,by="hh_size") %>% select(hh_size,bill,aff_grade,HBI,hh_income_upper_limit_Q1E,PPI) 


print(paste0("The Average Household Size for Naperville is ",
             round(pl.w$hh_size_avg_overallE,2)))

[1] “The Average Household Size for Naperville is 2.79”

print(paste0("At 50 GPCD, this corresponds to a monthly water volume of ",
             calculate_volume(hh_size=pl.w$hh_size_avg_overallE)))

[1] “At 50 GPCD, this corresponds to a monthly water volume of 4185”

print(paste0("This corresponds to an water and sewer bill of $",
             round(pl.w$bill_size_avg,2)))

[1] “This corresponds to an water and sewer bill of $66.67”

print(paste0("The Upper limit of the lowest quintile of annual household income for Naperville is $",
             pl.w$hh_income_upper_limit_Q1E,
             " or $",
             round(pl.w$hh_income_upper_limit_Q1E/12,2),
             " monthly"))

[1] “The Upper limit of the lowest quintile of annual household income for Naperville is $57958 or $4829.83 monthly”

print(paste0("The Poverty Prevalence Indicator (PPI) for Naperville overall is % ",
             round(pl.w$PPI,1)))

[1] “The Poverty Prevalence Indicator (PPI) for Naperville overall is % 9.4”

print(paste0("By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7"))

[1] “By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7”

datatable(pl.aff) %>%
  formatCurrency(c('bill','hh_income_upper_limit_Q1E'),'$') %>%
  formatRound(c('HBI','PPI'),2)

Below we visualize HBI, PPI, and the overall affordability grade as calcualted at the tract level.

# Tract
# Place
tr.w <- tr.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
tr.w <- tr.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
tr.w <- tr.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
tr.w <- tr.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
tr.w <- tr.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
tr.w <- tr.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
tr.w <- tr.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
tr.w <- tr.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))


tr.w$HBI_size_avg <- 100*tr.w$bill_size_avg/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_1 <- 100*tr.w$bill_size_1/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_2 <- 100*tr.w$bill_size_2/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_3 <- 100*tr.w$bill_size_3/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_4 <- 100*tr.w$bill_size_4/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_5 <- 100*tr.w$bill_size_5/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_6 <- 100*tr.w$bill_size_6/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_7 <- 100*tr.w$bill_size_7/(tr.w$hh_income_upper_limit_Q1E/12)


tr.w$PPI <- 100*(tr.w$inc_poverty_countE-tr.w$inc_poverty_count_gr2E)/tr.w$inc_poverty_countE
tr.w$aff_grade_avg <- gradeAffordability(tr.w$HBI_size_avg,tr.w$PPI)
tr.w$aff_grade_1 <- gradeAffordability(tr.w$HBI_size_1,tr.w$PPI)
tr.w$aff_grade_2 <- gradeAffordability(tr.w$HBI_size_2,tr.w$PPI)
tr.w$aff_grade_3 <- gradeAffordability(tr.w$HBI_size_3,tr.w$PPI)
tr.w$aff_grade_4 <- gradeAffordability(tr.w$HBI_size_4,tr.w$PPI)
tr.w$aff_grade_5 <- gradeAffordability(tr.w$HBI_size_5,tr.w$PPI)
tr.w$aff_grade_6 <- gradeAffordability(tr.w$HBI_size_6,tr.w$PPI)
tr.w$aff_grade_7 <- gradeAffordability(tr.w$HBI_size_7,tr.w$PPI)


hbi <- select(tr.w,
                  HBI_size_avg,
                  HBI_size_1,
                  HBI_size_2,
                  HBI_size_3,
                  HBI_size_4,
                  HBI_size_5,
                  HBI_size_6,
                  HBI_size_7,
                  )

ppi <- select(tr.w,
              PPI)

aff <- select(tr.w,
              aff_grade_avg,
              aff_grade_1,
              aff_grade_2,
              aff_grade_3,
              aff_grade_4,
              aff_grade_5,
              aff_grade_6,
              aff_grade_7)

m1 <- mapview(aff,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m2 <- mapview(hbi,burst=TRUE,col.regions=viridis::magma(7)) + mapview(b,alpha.regions=0,col.regions="green",stroke=TRUE,lwd=3,color="green",layer.name="Municipal Boundary")
m3 <- mapview(tr.w,zcol="hh_income_upper_limit_Q1E",col.regions=viridis::plasma(7),layer.name="Lower Quintile Income") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m4 <- mapview(ppi,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")


sync(m1,m2,m3,m4)

1.1.2, 1.1.3, 1.14, Cutoffs, Customer Assistance Programs (last 1 year), Conservation/ efficiency appliance incentive programs

Below, we simulate cutoff rates as a probabilistic function of HBI and PPI, and Customer Assistance participation rates as a probabilistic function of cutoff rates and white population share. We simulate participation in conservation/ efficiency appliance incentive programs as a function of homeownership rates and rate of population with at least a bachelors degree.

In reality, one would hope to use actually utility records of participation in these programs at the address level, or aggregated to census tracts.

tr.w$white_perc <- tr.w$pop_race_whiteE/tr.w$pop_race_countE

P1 <- select(tr.w, NAME, GEOID, aff_grade_avg, HBI_size_avg, PPI, hh_income_upper_limit_Q1E, hh_size_avg_overallE,white_perc, pop_owner_occupied_unitsE, pop_rental_unitsE, 
             pop_educ_attainment_countE, pop_educ_bachelorE, pop_educ_masterE, pop_educ_docE, pop_educ_profE) 

set.seed(20212002)

#Cutoffs and CAP
z <- 5  -0.03* P1$HBI_size_avg - 0.02*P1$PPI - 0.04 * P1$HBI_size_avg*P1$PPI
P1$Cutoff_Perc = 100*(1/(1+exp(z)))
P1$CAP_Perc <- P1$Cutoff_Perc * (1-0.3*(1-P1$white_perc))

#Conservation (ever)
P1$bachAbove <- 100*(P1$pop_educ_bachelorE + P1$pop_educ_masterE + P1$pop_educ_docE + P1$pop_educ_profE)/P1$pop_educ_attainment_countE
P1$own_rate <- 100*P1$pop_owner_occupied_unitsE/(P1$pop_rental_unitsE+P1$pop_owner_occupied_unitsE)

z <- 4 - 0.01 * P1$bachAbove - 0.02 * P1$own_rate
P1$Incentive_rate <- 100*(1/(1+exp(z)))

m5 <- mapview(P1,zcol="Cutoff_Perc", layer.name="Account cutoff (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m6 <- mapview(P1,zcol="CAP_Perc", layer.name="CAP participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m7 <- mapview(P1,zcol="Incentive_rate", layer.name="Efficiency Incentive participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")

sync(m5,m6,m7)

1.2 Water Quality (last 5 years)

We model all water quality events as basically random occurrences with probabilities between 0 and 10% and coverage of 50-100% at the tract level

1.2.1 Connections with Service Interruptions (not including nonpayment cutoffs or vacation-type shutoffs)

set.seed(32896)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$Interrup <- event*coverage

1.2.2 Connections with Boil water advisories

set.seed(328435396)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$boil_rate <- event*coverage

1.2.3 Health-based SDWA violations

This would be a Utility-Level measure only, and is available from EPA SDWIS

In this case, the API call for Naperville involves the PWSID IL0434670:

https://enviro.edap-cluster.com/cache/query/enviro3/efservice/VIOLATION/PWSID/IL0434670/IS_HEALTH_BASED_IND/Y/ROWS/0:100/JSON

In this case, there are no health-based violations in the past 10 years, so the value of the indicator is 0.

Pillar 2 | Economic and Workforce Development Contributions

Pillar 3: Resilience and Asset Management

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

  1. ↩︎